#include<vector>
#include<iostream>
#include<cmath>

using namespace std;
class Function
{
 public:
  double operator()(double _x);
};

class Newton_formula
{
 private:
  vector<double> x, y;
  vector<vector<double>> f;
 public:
  Newton_formula(
		 vector<double> _x,
		 vector<double> _y
		 )
    {
      x=_x;
      y=_y;
      f.resize(x.size()+1);
      for(int i =0;i<=x.size();i++)
	{
	  f[i].resize(x.size()+1);
	}
      for(int i=0;i < x.size();i++)
	{
	  f[i][0]=y[i];
	}
      for(int i =1; i < x.size();i++)
	{
	  for(int j= i; j< x.size();j++)
	    {
	      f[j][i]=(f[j][i-1]-f[j-1][i-1])/(x[j]-x[j-i]);
	    }
	}
	
    }
  vector<vector<double>> getf(){
    return f;
  }
  double calculate(double z )
  {
    double p=0;
    double k=1;
    for(int i=0; i<x.size();i++)
      {
	p = p + f[i][0] * k;
      k = 1;
      for(int j=0;j<=i;j++)
	  {
	    k=k*( z - x[j]);
	  }
      }
    return p;

  }
};


class Hermit
{
private:
  vector<double> x, y,x_;
  vector<vector<double>> f;
public:
  Hermit( vector<double> _x,
	  vector<double> _y,
	  vector<double> _x_){
    x=_x;
    y=_y;
    x_=_x_;
    f.resize(x.size()+1);
    for(int i =0;i<=x.size();i++)
	{
	  f[i].resize(x.size()+1);
	}
    for(int i=0;i < x.size();i++)
	{
	  f[i][0]=y[i];
	}
    for(int i=1;i < x.size();i = i + 2)
	{
	  f[i][1]=x_[i];
	}
    for(int i=2;i < x.size();i = i + 2)
	{
	  f[i][1]=(f[i][0]-f[i-1][0])/(x[i]-x[i-2]);
	}
    for(int i =2; i < x.size();i++)
	{
	  for(int j= i; j< x.size();j++)
	    {
	      f[j][i]=(f[j][i-1]-f[j-1][i-1])/(x[j]-x[j-i]);
	    }
	}
  }
  vector<vector<double>> getf(){
    return f;}
};
